from pylab import *
import scipy.optimize as sp
from pitubalib import *

#---------------------------------------------------------------------------
# PVT FROM SERGI RESERVOIR SAMPLE
#---------------------------------------------------------------------------

## PVT
## PVT do Bloco 3
fluido1=classPVT()
fluido1.pe=array([200.,170.,140.,110.,80.,60.,48.,35.,25.,15.,1.033])
fluido1.Boe=array([1.165,1.149,1.133,1.116,1.100,1.089,1.083,1.076,1.069,1.062,1.0000])
fluido1.Rse=array([66.9,57.8,48.6,39.4,30.2,24.1,20.4,16.4,13.1,9.4,0.0000])
fluido1.Bge=array([0.0047,0.055,0.0067,0.0085,0.0117,0.0156,0.0195,0.0310,0.0435,0.0720,1.0000])
fluido1.Bwe=1.0*ones(len(fluido1.pe))
fluido1.viscoe=array([5.81000,7.23000,8.05000,9.05000,10.23000,11.58000,12.17000,12.80000,14.12000,14.83000,15.57000,16.33000,17.12000,17.96000,18])
fluido1.viscwe=1.0*ones(len(fluido1.pe))
fluido1.co=93e-6
fluido1.Pb=48.
fluido1.T=273.+57.

## fluido simulado
fluido2=classPVT()
fluido2.pe=array([168.1700,140.1500,112.1200,84.0900,56.0600,50.4500,44.8500,39.2400,33.6300,28.0300,24.200,16.8200,11.2100,5.6100,1.00])
fluido2.Boe=array([1.1846,1.1616,1.1385,1.1155,1.0925,1.0874,1.0828,1.0786,1.0746,1.0705,1.0660,1.0608,1.0543,1.0457,1.0000])
fluido2.Rse=array([57.3700,48.0300,38.6867,29.3433,20.0000,18.1300,16.2300,14.3500,12.4700,10.5900,8.7200,6.8400,4.9200,2.9100,0.0000])
fluido2.Bge=array([0.0015,0.0022,0.0041,0.0076,0.0193,0.0216,0.0245,0.0283,0.0333,0.0404,0.0509,0.0686,0.1039,0.2097,1.0000])
fluido2.Bwe=1.03*ones(len(fluido2.pe))
fluido2.viscoe=array([5.81000,7.23000,8.05000,9.05000,10.23000,11.58000,12.17000,12.80000,14.12000,14.83000,15.57000,16.33000,17.12000,17.96000,18])
fluido2.viscwe=0.6*ones(len(fluido2.pe))
fluido2.co=1e-4
fluido2.Pb=35.0
fluido2.T=60.

fluido=fluido2

print gradient(fluido2.Boe)/gradient(fluido2.pe)


#---------------------------------------------------------------------------
# COMPUTING EQUATION COEFFICIENTS
# Np*(Bo2+(Rp-Rso2)*Bg2)+Wp*Bw2=N*(Eo+m*Eg+Efw)+(Wi+We)*Bw2+Gi*Bg2
#---------------------------------------------------------------------------
Eq=classEBM()
prod=genfromtxt('prod_AG_BL4.txt',delimiter='\t',skiprows=1)
Eq.Np=prod[:,0]*1e3
Eq.Gp=0.1*prod[:,1]*1e3
Eq.Wp=prod[:,2]*1e3
Eq.Wi=prod[:,3]*1e3
Eq.dt=prod[:,4]
Eq.We=Eq.Wp-Eq.Wi
Eq.Gi=0.*prod[:,1]

for i in range(1,len(Eq.We)):
	Eq.We[i]=max(Eq.We[i],Eq.We[i-1])

Eq.N=0.75e6

AG=classRes()
# DEFINING DATUM=-790M (AVERAGE RESERVOIR DEPTH)
# Pasta do poco CER-4: pe=91.3 kgf/cm2 @ MD=940m, COTA=-835.4,
AG.p0= 91. #91.3+0.08*(-835.4+800.)
AG.Pb0=fluido.Pb

AG.phi=0.21
AG.Sw0=0.13
AG.Swi=0.13
AG.Sor=0.25
AG.Vp0=Eq.N*fluido.Bo(AG.p0,AG.Pb0)/(1.-AG.Sw0-AG.Sg0)


# Petrophysics
AG.m=0.
AG.cr=2.5e-5 #3.5e-6 * 14.7/1.033#5e-5
petroPhys=classPetrophys()
Swe=array([0.13,0.23,0.34,0.44,0.54,0.65,0.75])
petroPhys.Swe=Swe
petroPhys.kroe=0.7*((1.-AG.Sor-Swe)/(1.-AG.Sw0-AG.Sor))**3
petroPhys.krwe=0.075*((Swe-AG.Sw0)/(1.-AG.Sw0-AG.Sor))**6

# Aquifer
aquifero=classAquifero()
aquifero.model='user_Np'
aquifero.Npe=Eq.Np
aquifero.Wee=Eq.We
aquifero.WW=15*Eq.N # 400e6*0.27 # 35e6
aquifero.p0=AG.p0
aquifero.cr=AG.cr
aquifero.cw=5.7e-5
aquifero.pvt=fluido
aquifero.k=250.
aquifero.L=1000.

print aquifero.We(20.,365,Np=100e3)

# TESTES DE PRESSAO

# DATUM=-800m
AG.teste.p=array([72.69,53.12,47.73,47.4,28.71,24.52,21.29,11.55,11.08,16.61,16.28,20.5])*1.
AG.teste.Np=array([3.87430,24.67000,36.40000,42.44000,65.32000,95.85000,108.40000,129.38000,144.10000,164.35000,171.40000,194.5])*1e3
#AG.teste.Gp=interp(AG.teste.Np,Eq.Np,Eq.Gp)
AG.teste.Wp=interp(AG.teste.Np,Eq.Np,Eq.Wp)
AG.teste.Wi=interp(AG.teste.Np,Eq.Np,Eq.Wi)

Eq.res=AG
Eq.aquifero=aquifero
Eq.pvt=fluido
Eq.petroPhys=petroPhys


#-----------------------------------------------------------------------
# RESULTADOS DO HISTORICO
#-----------------------------------------------------------------------
# COMPUTE GAS PRODUCTION FROM EXP PRESSURE
Gp=Eq.calcGp()
Eq.Gp=interp(Eq.Np,Gp[:,0],Gp[:,1])
savetxt('Gp.dat',Gp)

# COMPUTE MATERIAL BALANCE
out=classHist()
out=Eq.runHist(1.,AG.p0)
#Eq.Gp=zeros(len(Eq.Np))
#out=Eq.runHistcalcRGO(1.,AG.p0)

figure()
plot(Eq.Np*1e-6,out.p,'bo-')
plot(Eq.Np*1e-6,out.Pb,'--k')
errorbar(AG.teste.Np*1e-6,AG.teste.p,fmt='ro',yerr=60./14.23)
title('Material Balance :: CER Oil Field, Block: 4, Zone: AG.')
legend(('Calc Press','Bubble Press','Measured'))
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Reservoir Pressure (kgf/cm2)')
ylim([0,100])

figure()
plot(Eq.Np*1e-6,gradient(Eq.Gp)/gradient(Eq.Np),'rs-')
plot(Eq.Np*1e-6,zeros(len(Eq.Np)),'k:')
xlabel('Np (1e6 std m3)')
ylabel('GOR')
grid(1)

figure()
plot(Eq.Np*1e-6,Eq.Gp*1e-6,'bo-')
title('Material Balance :: CER Oil Field, Block: 4, Zone: Agua Grande.')
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Cumulative Gas production (1e6 std m3)')

figure()
plot(Eq.Np*1e-6,array(Eq.Wi)*1e-6,'bx-')
plot(Eq.Np*1e-6,array(out.We)*1e-6,'gx-')
title('Material Balance :: CER Oil Field, Block: 4, Zone: Agua Grande.')
legend(('Wi','We'),'upper left')
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Water Volume (1e6 std m3)')

figure()
plot(Eq.Np*1e-6,out.Sw,'bo-')
plot(Eq.Np*1e-6,out.So,'ko-')
plot(Eq.Np*1e-6,out.Sg,'ro-')
plot(Eq.Np*1e-6,array(out.Sg)+array(out.So)+array(out.Sw),'k-')
xlabel('Cumulative oil production (1e6 std m3)')
legend(('Sw','So','Sg','Sw+So+Sg'),'center left')
grid(1)
title('Saturations')

#figure()
#semilogy(range(len(Eq.Np)),abs(array(out.residual)),'k-') # Normalization due to sweep efficiency
#xlabel('step')
#ylabel('Residual')
#grid(1)

figure()
plot(1.-array(out.So),out.fw,'ro') # Normalization due to sweep efficiency
xlabel('1-So')
ylabel('BSW')
#xlim([0.,1.])
grid(1)

# Computing reservoir IP
#PI=gradient(Eq.Wp+Eq.Np)/array(out.p)*1e-3
#figure()
#plot(Eq.Np*1e-6,PI,'k-')
#xlabel('Cumulative oil production (1e6 std m3)')
#legend(('Reservoir PI'),'center left')
#grid(1)
#title('Reservoir PI')


#-----------------------------------------------------------------------
show()
print 'Gp max = ',(Eq.N*(fluido.Rs(AG.Pb0,AG.Pb0)-fluido.Rs(20.,AG.Pb0))+Eq.Np[-1]*fluido.Rs(20.,AG.Pb0))/1e6
print 'Porous Volume = ',AG.Vp0/1e6
print 'Savg = ',AG.Vp0/18.
print 'hgas = ',AG.Vp0*out.Sg[-1]/18.
